Coral Wound Healing Analysis: Effects of Damselfish Presence and Wound Size

Overview

Study Design

Research Question: Do damselfish (Dascyllus flavicaudus) accelerate coral wound healing rates in Pocillopora spp.?

Experimental Design:

• Factors: Fish presence (Present/Absent) × Wound size (None/Small/Large)
• Sample size: 72 coral fragments
• Duration: 21 days
• Location: Mo’orea, French Polynesia
• Response variables:
  • Healing rate (cm²/day)
  • Percent wound closure
  • Binary healing outcome

Analysis Overview

This document contains:

1. Data Import & Cleaning
2. Exploratory Data Analysis
3. Mixed-Effects Models (healing rate)
4. Publication Figure Generation
5. Model Diagnostics

---

1. Setup & Data Import

Load Required Packages

# ==================== Core Packages ====================
library(tidyverse)      # Data wrangling and visualization
library(here)           # Relative file paths
library(janitor)        # Data cleaning

# ==================== Statistical Modeling ====================
library(lme4)           # Linear mixed-effects models
library(glmmTMB)        # Generalized linear mixed models
library(DHARMa)         # Residual diagnostics for GLMMs
library(MuMIn)          # Model selection
library(parameters)     # Model parameter extraction
library(broom.mixed)    # Tidy model outputs
library(ggeffects)      # Marginal effects

# ==================== Visualization ====================
library(ggpubr)         # Publication-ready plots
library(viridis)        # Color palettes
library(patchwork)      # Combine plots
library(gt)             # Publication-quality tables

# ==================== Fonts ====================
library(sysfonts)       # Font management
library(showtext)       # Font rendering

# ==================== Optional: Advanced Models ====================
if (!requireNamespace("gamlss", quietly = TRUE)) {
  install.packages("gamlss")
}
library(gamlss)
library(gamlss.dist)

# Setup font rendering
font_add_google(name = "Josefin Sans", family = "josefin")
showtext_auto()

# Set default theme
theme_set(theme_bw() + theme(panel.grid = element_blank(), plot.background = element_blank()))

Define Color Palettes

# Consistent color schemes for all figures
wound_palette <- c(
  "Large" = "#FF8B00",     # Orange
  "Small" = "#be8333",     # Tan
  "No Wound" = "#485900"   # Olive green
)

fish_palette <- c(
  "Fish" = "#2b4b52",      # Teal
  "No Fish" = "#ad301a"    # Burnt red
)

Import Raw Data

# Import wound closure binary data (healed yes/no)
wound_closure <- read_csv(
  here("data", "fish_regen_mastersheet_wound closure_necrosis_sa - mastersheet.csv")
) %>%
  clean_names()

# Import quantitative healing rate data
percent_closure <- read_csv(
  here("data", "final_wound_size_dascyllus_project - Sheet1.csv")
) %>%
  clean_names()

cat("✓ Data imported successfully\n")

## ✓ Data imported successfully

cat("  - Wound closure binary data:", nrow(wound_closure), "observations\n")

##   - Wound closure binary data: 72 observations

cat("  - Percent closure data:", nrow(percent_closure), "observations\n")

##   - Percent closure data: 48 observations

---

2. Data Preparation

Treatment Coding Function

#' Convert numeric treatment codes to labeled factors
#'
#' @param df Data frame containing fish and wound columns
#' @return Data frame with properly ordered factor levels
prep_treatments <- function(df) {
  df %>%
    mutate(
      # Convert fish coding: 1 = Fish present, 0 = Fish absent
      fish = ifelse(fish == 1, "Fish", ifelse(fish == 0, "No Fish", fish)),

      # Convert wound size coding: 0 = None, 1 = Small, 2 = Large
      wound = case_when(
        wound == 0 ~ "No Wound",
        wound == 1 ~ "Small",
        wound == 2 ~ "Large",
        TRUE ~ as.character(wound)
      )
    ) %>%
    mutate(
      # Set factor levels in logical order
      fish  = factor(fish, levels = c("No Fish", "Fish")),
      wound = factor(wound, levels = c("No Wound", "Small", "Large"))
    )
}

Clean Wound Closure Data

# Process binary wound closure data
wound_closure_clean <- wound_closure %>%
  filter(!is.na(wound_close), wound_close != "") %>%
  # Standardize "barely healed" responses to "no"
  mutate(wound_close = ifelse(wound_close %in% c("no-barely", "no - barely"), "no", wound_close)) %>%
  mutate(wound_close = ifelse(wound_close == "yes", "Yes", "No")) %>%
  prep_treatments() %>%
  mutate(wound_close = as.character(wound_close))

# Create binary version for GLM modeling
wound_closure_glm <- wound_closure %>%
  filter(!is.na(wound_close), wound_close != "") %>%
  mutate(wound_bin = ifelse(wound_close == "yes", 1, 0)) %>%
  prep_treatments() %>%
  filter(!is.na(fish), !is.na(wound), !is.na(tank)) %>%
  droplevels()

# Quality check: Ensure factors have adequate levels
if (nlevels(wound_closure_glm$fish) < 2 | nlevels(wound_closure_glm$wound) < 2) {
  stop("ERROR: Fish or wound factors have <2 levels. Check input data.")
}

cat("✓ Data cleaning complete\n")

## ✓ Data cleaning complete

cat("  - Binary healing data:", nrow(wound_closure_clean), "corals\n")

##   - Binary healing data: 48 corals

Clean Healing Rate Data

# Process continuous healing rate data
percent_closure <- percent_closure %>%
  filter(percent_healed >= 0) %>%  # Remove invalid percentages
  prep_treatments() %>%
  mutate(tank = as.factor(tank))

# Convert percent to proportion (0-1 scale)
percent_closure$prop <- percent_closure$percent_healed / 100
percent_closure$prop_remaining <- 1 - percent_closure$prop

cat("✓ Healing rate data prepared\n")

## ✓ Healing rate data prepared

cat("  - Valid observations:", nrow(percent_closure), "corals\n")

##   - Valid observations: 47 corals

cat("  - Mean healing rate:", round(mean(percent_closure$healing_rate, na.rm = TRUE), 3), "cm²/day\n")

##   - Mean healing rate: 0.119 cm²/day

Data Summary

# Cross-tabulation of treatments
cat("\n=== Treatment Sample Sizes ===\n")

## 
## === Treatment Sample Sizes ===

print(table(percent_closure$fish, percent_closure$wound))

##          
##           No Wound Small Large
##   No Fish        0    11    12
##   Fish           0    12    12

# Calculate healing proportions by treatment
wound_prop <- wound_closure_clean %>%
  group_by(wound, fish) %>%
  summarize(
    n = n(),
    prop_healed = mean(wound_close == "Yes"),
    .groups = "drop"
  )

cat("\n=== Proportion Healed by Treatment ===\n")

## 
## === Proportion Healed by Treatment ===

print(wound_prop)

## # A tibble: 4 × 4
##   wound fish        n prop_healed
##   <fct> <fct>   <int>       <dbl>
## 1 Small No Fish    12       0.583
## 2 Small Fish       12       0.917
## 3 Large No Fish    12       0.417
## 4 Large Fish       12       0.75

---

3. Wound Size Characterization

# Summarize initial wound sizes
wounds_avgs <- percent_closure %>%
  group_by(wound) %>%
  summarize(
    mean_size_cm2 = mean(initial_wound_size_cm, na.rm = TRUE),
    sd_size_cm2 = sd(initial_wound_size_cm, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  ) %>%
  filter(wound != "No Wound")

# Display as formatted table
wounds_avgs %>%
  gt() %>%
  tab_header(title = "Initial Wound Size by Treatment") %>%
  fmt_number(columns = c(mean_size_cm2, sd_size_cm2), decimals = 2) %>%
  cols_label(
    wound = "Wound Size",
    mean_size_cm2 = "Mean (cm²)",
    sd_size_cm2 = "SD (cm²)",
    n = "N"
  )

Initial Wound Size by Treatment
Wound Size	Mean (cm²)	SD (cm²)	N
Small	1.48	0.40	23
Large	3.82	1.08	24

---

4. Statistical Models: Healing Rate

Model Description

We fit linear mixed-effects models to test how fish presence and wound size affect healing rate (continuous response variable in cm²/day).

Model structure: - Fixed effects: Fish presence, Wound size, and their interaction - Random effect: Tank (accounts for shared tank environment) - Response: Healing rate (cm²/day)

Fit Candidate Models

# Full model with interaction
rate_full <- lmer(
  healing_rate ~ fish * wound + (1 | tank),
  data = percent_closure
)

# Additive model (no interaction)
rate_no_interaction <- lmer(
  healing_rate ~ fish + wound + (1 | tank),
  data = percent_closure
)

# Fish-only model
rate_no_fish <- lmer(
  healing_rate ~ wound + (1 | tank),
  data = percent_closure
)

# Wound-only model
rate_no_wound <- lmer(
  healing_rate ~ fish + (1 | tank),
  data = percent_closure
)

cat("✓ Models fitted successfully\n")

## ✓ Models fitted successfully

Model Comparison: Likelihood Ratio Tests

# Test for interaction effect
lrt_interaction_zero_rate <- anova(rate_no_interaction, rate_full, test = "Chisq")

# Test for fish main effect
lrt_fish_zero_rate <- anova(rate_no_fish, rate_no_interaction, test = "Chisq")

# Test for wound main effect
lrt_wound_zero_rate <- anova(rate_no_wound, rate_no_interaction, test = "Chisq")

# Compile results into table
lrt_table_rate <- tibble(
  Test = c("Interaction (Fish × Wound)", "Fish Effect", "Wound Effect"),
  ChiSq = c(
    lrt_interaction_zero_rate$Chisq[2],
    lrt_fish_zero_rate$Chisq[2],
    lrt_wound_zero_rate$Chisq[2]
  ),
  Df = c(
    lrt_interaction_zero_rate$Df[2],
    lrt_fish_zero_rate$Df[2],
    lrt_wound_zero_rate$Df[2]
  ),
  p_value = c(
    lrt_interaction_zero_rate$`Pr(>Chisq)`[2],
    lrt_fish_zero_rate$`Pr(>Chisq)`[2],
    lrt_wound_zero_rate$`Pr(>Chisq)`[2]
  )
) %>%
  mutate(across(where(is.numeric), round, 3))

# Format as publication table
gt_table_rate <- lrt_table_rate %>%
  gt() %>%
  tab_header(
    title = "Likelihood Ratio Tests for Effects on Healing Rate"
  ) %>%
  cols_label(
    Test = "Model Comparison",
    ChiSq = "χ²",
    Df = "df",
    p_value = "P-value"
  ) %>%
  fmt_number(columns = c(ChiSq, p_value), decimals = 3) %>%
  tab_options(
    table.font.size = 14,
    heading.title.font.size = 16
  )

# Display table
print(gt_table_rate)

## <div id="ldbicnuyvi" style="padding-left:0px;padding-right:0px;padding-top:10px;padding-bottom:10px;overflow-x:auto;overflow-y:auto;width:auto;height:auto;">
##   <style>#ldbicnuyvi table {
##   font-family: system-ui, 'Segoe UI', Roboto, Helvetica, Arial, sans-serif, 'Apple Color Emoji', 'Segoe UI Emoji', 'Segoe UI Symbol', 'Noto Color Emoji';
##   -webkit-font-smoothing: antialiased;
##   -moz-osx-font-smoothing: grayscale;
## }
## 
## #ldbicnuyvi thead, #ldbicnuyvi tbody, #ldbicnuyvi tfoot, #ldbicnuyvi tr, #ldbicnuyvi td, #ldbicnuyvi th {
##   border-style: none;
## }
## 
## #ldbicnuyvi p {
##   margin: 0;
##   padding: 0;
## }
## 
## #ldbicnuyvi .gt_table {
##   display: table;
##   border-collapse: collapse;
##   line-height: normal;
##   margin-left: auto;
##   margin-right: auto;
##   color: #333333;
##   font-size: 14px;
##   font-weight: normal;
##   font-style: normal;
##   background-color: #FFFFFF;
##   width: auto;
##   border-top-style: solid;
##   border-top-width: 2px;
##   border-top-color: #A8A8A8;
##   border-right-style: none;
##   border-right-width: 2px;
##   border-right-color: #D3D3D3;
##   border-bottom-style: solid;
##   border-bottom-width: 2px;
##   border-bottom-color: #A8A8A8;
##   border-left-style: none;
##   border-left-width: 2px;
##   border-left-color: #D3D3D3;
## }
## 
## #ldbicnuyvi .gt_caption {
##   padding-top: 4px;
##   padding-bottom: 4px;
## }
## 
## #ldbicnuyvi .gt_title {
##   color: #333333;
##   font-size: 16px;
##   font-weight: initial;
##   padding-top: 4px;
##   padding-bottom: 4px;
##   padding-left: 5px;
##   padding-right: 5px;
##   border-bottom-color: #FFFFFF;
##   border-bottom-width: 0;
## }
## 
## #ldbicnuyvi .gt_subtitle {
##   color: #333333;
##   font-size: 85%;
##   font-weight: initial;
##   padding-top: 3px;
##   padding-bottom: 5px;
##   padding-left: 5px;
##   padding-right: 5px;
##   border-top-color: #FFFFFF;
##   border-top-width: 0;
## }
## 
## #ldbicnuyvi .gt_heading {
##   background-color: #FFFFFF;
##   text-align: center;
##   border-bottom-color: #FFFFFF;
##   border-left-style: none;
##   border-left-width: 1px;
##   border-left-color: #D3D3D3;
##   border-right-style: none;
##   border-right-width: 1px;
##   border-right-color: #D3D3D3;
## }
## 
## #ldbicnuyvi .gt_bottom_border {
##   border-bottom-style: solid;
##   border-bottom-width: 2px;
##   border-bottom-color: #D3D3D3;
## }
## 
## #ldbicnuyvi .gt_col_headings {
##   border-top-style: solid;
##   border-top-width: 2px;
##   border-top-color: #D3D3D3;
##   border-bottom-style: solid;
##   border-bottom-width: 2px;
##   border-bottom-color: #D3D3D3;
##   border-left-style: none;
##   border-left-width: 1px;
##   border-left-color: #D3D3D3;
##   border-right-style: none;
##   border-right-width: 1px;
##   border-right-color: #D3D3D3;
## }
## 
## #ldbicnuyvi .gt_col_heading {
##   color: #333333;
##   background-color: #FFFFFF;
##   font-size: 100%;
##   font-weight: normal;
##   text-transform: inherit;
##   border-left-style: none;
##   border-left-width: 1px;
##   border-left-color: #D3D3D3;
##   border-right-style: none;
##   border-right-width: 1px;
##   border-right-color: #D3D3D3;
##   vertical-align: bottom;
##   padding-top: 5px;
##   padding-bottom: 6px;
##   padding-left: 5px;
##   padding-right: 5px;
##   overflow-x: hidden;
## }
## 
## #ldbicnuyvi .gt_column_spanner_outer {
##   color: #333333;
##   background-color: #FFFFFF;
##   font-size: 100%;
##   font-weight: normal;
##   text-transform: inherit;
##   padding-top: 0;
##   padding-bottom: 0;
##   padding-left: 4px;
##   padding-right: 4px;
## }
## 
## #ldbicnuyvi .gt_column_spanner_outer:first-child {
##   padding-left: 0;
## }
## 
## #ldbicnuyvi .gt_column_spanner_outer:last-child {
##   padding-right: 0;
## }
## 
## #ldbicnuyvi .gt_column_spanner {
##   border-bottom-style: solid;
##   border-bottom-width: 2px;
##   border-bottom-color: #D3D3D3;
##   vertical-align: bottom;
##   padding-top: 5px;
##   padding-bottom: 5px;
##   overflow-x: hidden;
##   display: inline-block;
##   width: 100%;
## }
## 
## #ldbicnuyvi .gt_spanner_row {
##   border-bottom-style: hidden;
## }
## 
## #ldbicnuyvi .gt_group_heading {
##   padding-top: 8px;
##   padding-bottom: 8px;
##   padding-left: 5px;
##   padding-right: 5px;
##   color: #333333;
##   background-color: #FFFFFF;
##   font-size: 100%;
##   font-weight: initial;
##   text-transform: inherit;
##   border-top-style: solid;
##   border-top-width: 2px;
##   border-top-color: #D3D3D3;
##   border-bottom-style: solid;
##   border-bottom-width: 2px;
##   border-bottom-color: #D3D3D3;
##   border-left-style: none;
##   border-left-width: 1px;
##   border-left-color: #D3D3D3;
##   border-right-style: none;
##   border-right-width: 1px;
##   border-right-color: #D3D3D3;
##   vertical-align: middle;
##   text-align: left;
## }
## 
## #ldbicnuyvi .gt_empty_group_heading {
##   padding: 0.5px;
##   color: #333333;
##   background-color: #FFFFFF;
##   font-size: 100%;
##   font-weight: initial;
##   border-top-style: solid;
##   border-top-width: 2px;
##   border-top-color: #D3D3D3;
##   border-bottom-style: solid;
##   border-bottom-width: 2px;
##   border-bottom-color: #D3D3D3;
##   vertical-align: middle;
## }
## 
## #ldbicnuyvi .gt_from_md > :first-child {
##   margin-top: 0;
## }
## 
## #ldbicnuyvi .gt_from_md > :last-child {
##   margin-bottom: 0;
## }
## 
## #ldbicnuyvi .gt_row {
##   padding-top: 8px;
##   padding-bottom: 8px;
##   padding-left: 5px;
##   padding-right: 5px;
##   margin: 10px;
##   border-top-style: solid;
##   border-top-width: 1px;
##   border-top-color: #D3D3D3;
##   border-left-style: none;
##   border-left-width: 1px;
##   border-left-color: #D3D3D3;
##   border-right-style: none;
##   border-right-width: 1px;
##   border-right-color: #D3D3D3;
##   vertical-align: middle;
##   overflow-x: hidden;
## }
## 
## #ldbicnuyvi .gt_stub {
##   color: #333333;
##   background-color: #FFFFFF;
##   font-size: 100%;
##   font-weight: initial;
##   text-transform: inherit;
##   border-right-style: solid;
##   border-right-width: 2px;
##   border-right-color: #D3D3D3;
##   padding-left: 5px;
##   padding-right: 5px;
## }
## 
## #ldbicnuyvi .gt_stub_row_group {
##   color: #333333;
##   background-color: #FFFFFF;
##   font-size: 100%;
##   font-weight: initial;
##   text-transform: inherit;
##   border-right-style: solid;
##   border-right-width: 2px;
##   border-right-color: #D3D3D3;
##   padding-left: 5px;
##   padding-right: 5px;
##   vertical-align: top;
## }
## 
## #ldbicnuyvi .gt_row_group_first td {
##   border-top-width: 2px;
## }
## 
## #ldbicnuyvi .gt_row_group_first th {
##   border-top-width: 2px;
## }
## 
## #ldbicnuyvi .gt_summary_row {
##   color: #333333;
##   background-color: #FFFFFF;
##   text-transform: inherit;
##   padding-top: 8px;
##   padding-bottom: 8px;
##   padding-left: 5px;
##   padding-right: 5px;
## }
## 
## #ldbicnuyvi .gt_first_summary_row {
##   border-top-style: solid;
##   border-top-color: #D3D3D3;
## }
## 
## #ldbicnuyvi .gt_first_summary_row.thick {
##   border-top-width: 2px;
## }
## 
## #ldbicnuyvi .gt_last_summary_row {
##   padding-top: 8px;
##   padding-bottom: 8px;
##   padding-left: 5px;
##   padding-right: 5px;
##   border-bottom-style: solid;
##   border-bottom-width: 2px;
##   border-bottom-color: #D3D3D3;
## }
## 
## #ldbicnuyvi .gt_grand_summary_row {
##   color: #333333;
##   background-color: #FFFFFF;
##   text-transform: inherit;
##   padding-top: 8px;
##   padding-bottom: 8px;
##   padding-left: 5px;
##   padding-right: 5px;
## }
## 
## #ldbicnuyvi .gt_first_grand_summary_row {
##   padding-top: 8px;
##   padding-bottom: 8px;
##   padding-left: 5px;
##   padding-right: 5px;
##   border-top-style: double;
##   border-top-width: 6px;
##   border-top-color: #D3D3D3;
## }
## 
## #ldbicnuyvi .gt_last_grand_summary_row_top {
##   padding-top: 8px;
##   padding-bottom: 8px;
##   padding-left: 5px;
##   padding-right: 5px;
##   border-bottom-style: double;
##   border-bottom-width: 6px;
##   border-bottom-color: #D3D3D3;
## }
## 
## #ldbicnuyvi .gt_striped {
##   background-color: rgba(128, 128, 128, 0.05);
## }
## 
## #ldbicnuyvi .gt_table_body {
##   border-top-style: solid;
##   border-top-width: 2px;
##   border-top-color: #D3D3D3;
##   border-bottom-style: solid;
##   border-bottom-width: 2px;
##   border-bottom-color: #D3D3D3;
## }
## 
## #ldbicnuyvi .gt_footnotes {
##   color: #333333;
##   background-color: #FFFFFF;
##   border-bottom-style: none;
##   border-bottom-width: 2px;
##   border-bottom-color: #D3D3D3;
##   border-left-style: none;
##   border-left-width: 2px;
##   border-left-color: #D3D3D3;
##   border-right-style: none;
##   border-right-width: 2px;
##   border-right-color: #D3D3D3;
## }
## 
## #ldbicnuyvi .gt_footnote {
##   margin: 0px;
##   font-size: 90%;
##   padding-top: 4px;
##   padding-bottom: 4px;
##   padding-left: 5px;
##   padding-right: 5px;
## }
## 
## #ldbicnuyvi .gt_sourcenotes {
##   color: #333333;
##   background-color: #FFFFFF;
##   border-bottom-style: none;
##   border-bottom-width: 2px;
##   border-bottom-color: #D3D3D3;
##   border-left-style: none;
##   border-left-width: 2px;
##   border-left-color: #D3D3D3;
##   border-right-style: none;
##   border-right-width: 2px;
##   border-right-color: #D3D3D3;
## }
## 
## #ldbicnuyvi .gt_sourcenote {
##   font-size: 90%;
##   padding-top: 4px;
##   padding-bottom: 4px;
##   padding-left: 5px;
##   padding-right: 5px;
## }
## 
## #ldbicnuyvi .gt_left {
##   text-align: left;
## }
## 
## #ldbicnuyvi .gt_center {
##   text-align: center;
## }
## 
## #ldbicnuyvi .gt_right {
##   text-align: right;
##   font-variant-numeric: tabular-nums;
## }
## 
## #ldbicnuyvi .gt_font_normal {
##   font-weight: normal;
## }
## 
## #ldbicnuyvi .gt_font_bold {
##   font-weight: bold;
## }
## 
## #ldbicnuyvi .gt_font_italic {
##   font-style: italic;
## }
## 
## #ldbicnuyvi .gt_super {
##   font-size: 65%;
## }
## 
## #ldbicnuyvi .gt_footnote_marks {
##   font-size: 75%;
##   vertical-align: 0.4em;
##   position: initial;
## }
## 
## #ldbicnuyvi .gt_asterisk {
##   font-size: 100%;
##   vertical-align: 0;
## }
## 
## #ldbicnuyvi .gt_indent_1 {
##   text-indent: 5px;
## }
## 
## #ldbicnuyvi .gt_indent_2 {
##   text-indent: 10px;
## }
## 
## #ldbicnuyvi .gt_indent_3 {
##   text-indent: 15px;
## }
## 
## #ldbicnuyvi .gt_indent_4 {
##   text-indent: 20px;
## }
## 
## #ldbicnuyvi .gt_indent_5 {
##   text-indent: 25px;
## }
## 
## #ldbicnuyvi .katex-display {
##   display: inline-flex !important;
##   margin-bottom: 0.75em !important;
## }
## 
## #ldbicnuyvi div.Reactable > div.rt-table > div.rt-thead > div.rt-tr.rt-tr-group-header > div.rt-th-group:after {
##   height: 0px !important;
## }
## </style>
##   <table class="gt_table" data-quarto-disable-processing="false" data-quarto-bootstrap="false">
##   <thead>
##     <tr class="gt_heading">
##       <td colspan="4" class="gt_heading gt_title gt_font_normal gt_bottom_border" style>Likelihood Ratio Tests for Effects on Healing Rate</td>
##     </tr>
##     
##     <tr class="gt_col_headings">
##       <th class="gt_col_heading gt_columns_bottom_border gt_left" rowspan="1" colspan="1" scope="col" id="Test">Model Comparison</th>
##       <th class="gt_col_heading gt_columns_bottom_border gt_right" rowspan="1" colspan="1" scope="col" id="ChiSq">χ²</th>
##       <th class="gt_col_heading gt_columns_bottom_border gt_right" rowspan="1" colspan="1" scope="col" id="Df">df</th>
##       <th class="gt_col_heading gt_columns_bottom_border gt_right" rowspan="1" colspan="1" scope="col" id="p_value">P-value</th>
##     </tr>
##   </thead>
##   <tbody class="gt_table_body">
##     <tr><td headers="Test" class="gt_row gt_left">Interaction (Fish × Wound)</td>
## <td headers="ChiSq" class="gt_row gt_right">4.416</td>
## <td headers="Df" class="gt_row gt_right">1</td>
## <td headers="p_value" class="gt_row gt_right">0.036</td></tr>
##     <tr><td headers="Test" class="gt_row gt_left">Fish Effect</td>
## <td headers="ChiSq" class="gt_row gt_right">8.047</td>
## <td headers="Df" class="gt_row gt_right">1</td>
## <td headers="p_value" class="gt_row gt_right">0.005</td></tr>
##     <tr><td headers="Test" class="gt_row gt_left">Wound Effect</td>
## <td headers="ChiSq" class="gt_row gt_right">40.900</td>
## <td headers="Df" class="gt_row gt_right">1</td>
## <td headers="p_value" class="gt_row gt_right">0.000</td></tr>
##   </tbody>
##   
##   
## </table>
## </div>

# Save table
gtsave(gt_table_rate, here("figures", "lrt_table_rate.html"))

Model Summary

# Display full model summary
summary(rate_full)

## Linear mixed model fit by REML ['lmerMod']
## Formula: healing_rate ~ fish * wound + (1 | tank)
##    Data: percent_closure
## 
## REML criterion at convergence: -140
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.72955 -0.53791  0.04203  0.56225  2.13584 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  tank     (Intercept) 6.678e-05 0.008172
##  Residual             1.752e-03 0.041852
## Number of obs: 47, groups:  tank, 6
## 
## Fixed effects:
##                     Estimate Std. Error t value
## (Intercept)          0.05682    0.01348   4.214
## fishFish             0.02166    0.01871   1.158
## woundLarge           0.07441    0.01748   4.257
## fishFish:woundLarge  0.05077    0.02444   2.077
## 
## Correlation of Fixed Effects:
##             (Intr) fshFsh wndLrg
## fishFish    -0.721              
## woundLarge  -0.677  0.488       
## fshFsh:wndL  0.484 -0.668 -0.715

# Extract and display fixed effects
cat("\n=== Fixed Effects (Full Model) ===\n")

## 
## === Fixed Effects (Full Model) ===

print(tidy(rate_full, effects = "fixed", conf.int = TRUE))

## # A tibble: 4 × 7
##   effect term                estimate std.error statistic conf.low conf.high
##   <chr>  <chr>                  <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 fixed  (Intercept)           0.0568    0.0135      4.21  0.0304     0.0832
## 2 fixed  fishFish              0.0217    0.0187      1.16 -0.0150     0.0583
## 3 fixed  woundLarge            0.0744    0.0175      4.26  0.0402     0.109 
## 4 fixed  fishFish:woundLarge   0.0508    0.0244      2.08  0.00287    0.0987

---

5. Publication Figure: Healing Rate

# Create publication-quality figure
healing_plot <- percent_closure %>%
  ggplot(aes(x = wound, y = healing_rate, color = fish, shape = fish)) +

  # Individual data points (jittered and dodged)
  geom_jitter(
    aes(color = fish, shape = fish),
    position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.6),
    size = 2.5,
    alpha = 0.25,
    shape = 16
  ) +

  # Mean points
  stat_summary(
    fun = mean,
    geom = "point",
    size = 5,
    position = position_dodge(width = 0.6)
  ) +

  # Error bars (mean ± SD)
  stat_summary(
    fun.data = function(x) {
      data.frame(
        y = mean(x),
        ymin = mean(x) - sd(x),
        ymax = mean(x) + sd(x)
      )
    },
    geom = "errorbar",
    width = 0,
    size = 1,
    position = position_dodge(width = 0.6)
  ) +

  # Color and shape scales
  scale_color_manual(
    values = fish_palette,
    name = "Fish Presence",
    labels = c("No Fish", "Fish")
  ) +
  scale_shape_manual(
    values = c(16, 18),
    name = "Fish Presence",
    labels = c("No Fish", "Fish")
  )+
  scale_y_continuous(limits = c(0, 0.35)) +

  # Theme and aesthetics
  theme_minimal(base_size = 16) +
  theme(
    panel.grid = element_blank(),
    axis.title = element_text(size = 20, face = "bold"),
    axis.text.x = element_text(size = 18),
    axis.text.y = element_text(size = 16),
    legend.position = "top",
    legend.title = element_text(size = 18, face = "bold"),
    legend.text = element_text(size = 16),
    legend.key.size = unit(1.5, "lines"),
    legend.margin = margin(0, 0, 10, 0),
    legend.box.spacing = unit(0.5, "lines"),
    axis.line.x = element_line(color = "black", size = 0.5),
    axis.line.y = element_line(color = "black", size = 0.5)
  ) +
  labs(
    y = "Healing Rate (cm²/day)",
    x = "Wound Size"
  )

# Display plot
print(healing_plot)

# Save as PDF (vector, print-quality)
ggsave(
  filename = here("figures", "wound_closure", "summary_line_healing_rate_by_wound_fish.pdf"),
  plot = healing_plot,
  width = 7,
  height = 5,
  dpi = 600,
  device = cairo_pdf
)

# Save as PNG (raster, screen-quality)
ggsave(
  filename = here("figures", "wound_closure", "summary_line_healing_rate_by_wound_fish.png"),
  plot = healing_plot,
  width = 8,
  height = 6,
  dpi = 300,
  bg = "white"
)

cat("✓ Figure saved to figures/wound_closure/\n")

## ✓ Figure saved to figures/wound_closure/

---

6. Summary Statistics by Treatment

# Wound size effect
rate_avg_wound <- percent_closure %>%
  group_by(wound) %>%
  summarize(
    mean_rate = mean(healing_rate, na.rm = TRUE),
    sd_rate = sd(healing_rate, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  )

cat("=== Healing Rate by Wound Size ===\n")

## === Healing Rate by Wound Size ===

print(rate_avg_wound)

## # A tibble: 2 × 4
##   wound mean_rate sd_rate     n
##   <fct>     <dbl>   <dbl> <int>
## 1 Small    0.0682  0.0199    23
## 2 Large    0.167   0.0669    24

# Fish presence effect
rate_avg_fish <- percent_closure %>%
  group_by(fish) %>%
  summarize(
    mean_rate = mean(healing_rate, na.rm = TRUE),
    sd_rate = sd(healing_rate, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  )

cat("\n=== Healing Rate by Fish Presence ===\n")

## 
## === Healing Rate by Fish Presence ===

print(rate_avg_fish)

## # A tibble: 2 × 4
##   fish    mean_rate sd_rate     n
##   <fct>       <dbl>   <dbl> <int>
## 1 No Fish    0.0957  0.0566    23
## 2 Fish       0.141   0.0759    24

# Interaction: Fish × Wound
rate_avg_interaction <- percent_closure %>%
  group_by(fish, wound) %>%
  summarize(
    mean_rate = mean(healing_rate, na.rm = TRUE),
    sd_rate = sd(healing_rate, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  )

cat("\n=== Healing Rate by Fish × Wound ===\n")

## 
## === Healing Rate by Fish × Wound ===

print(rate_avg_interaction)

## # A tibble: 4 × 5
##   fish    wound mean_rate sd_rate     n
##   <fct>   <fct>     <dbl>   <dbl> <int>
## 1 No Fish Small    0.0570  0.0114    11
## 2 No Fish Large    0.131   0.0586    12
## 3 Fish    Small    0.0785  0.0208    12
## 4 Fish    Large    0.204   0.0553    12

---

7. Model Diagnostics

DHARMa Residual Diagnostics

library(DHARMa)

# Simulate residuals for full model
sim_rate_full <- simulateResiduals(fittedModel = rate_full, n = 1000)

# Save DHARMa diagnostic plot
png(
  here("figures", "diagnostics", "sim_rate_full.png"),
  width = 7, height = 6, units = "in", res = 600, bg = "white"
)
plot(sim_rate_full, main = "DHARMa Diagnostics: Full Model")
dev.off()

## quartz_off_screen 
##                 2

# Statistical tests
cat("\n=== Dispersion Test ===\n")

## 
## === Dispersion Test ===

print(testDispersion(sim_rate_full))

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.92535, p-value = 0.754
## alternative hypothesis: two.sided

cat("\n=== Zero-Inflation Test ===\n")

## 
## === Zero-Inflation Test ===

print(testZeroInflation(sim_rate_full))

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided

Diagnostics for All Candidate Models

model_list <- list(
  Full = rate_full,
  NoInteraction = rate_no_interaction,
  NoFish = rate_no_fish,
  NoWound = rate_no_wound
)

for (nm in names(model_list)) {
  cat("\n========== Diagnostics for", nm, "model ==========\n")
  sim <- simulateResiduals(model_list[[nm]], n = 1000)

  # Save DHARMa plot
  png(
    here("figures", "diagnostics", paste0("resid_", nm, ".png")),
    width = 7, height = 6, units = "in", res = 600, bg = "white"
  )
  plot(sim, main = paste("DHARMa Diagnostics:", nm, "Model"))
  dev.off()

  # Print test results
  cat("\n  Dispersion test:\n")
  print(testDispersion(sim))
  cat("\n  Outlier test:\n")
  print(testOutliers(sim))
}

## 
## ========== Diagnostics for Full model ==========

## 
##   Dispersion test:

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.92535, p-value = 0.754
## alternative hypothesis: two.sided
## 
## 
##   Outlier test:

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  sim
## outliers at both margin(s) = 0, observations = 47, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.00000000 0.07548573
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                                    0 
## 
## 
## ========== Diagnostics for NoInteraction model ==========

## 
##   Dispersion test:

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.94733, p-value = 0.842
## alternative hypothesis: two.sided
## 
## 
##   Outlier test:

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  sim
## outliers at both margin(s) = 1, observations = 47, p-value = 0.08972
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0005385317 0.1129377171
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                            0.0212766 
## 
## 
## ========== Diagnostics for NoFish model ==========

## 
##   Dispersion test:

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.97263, p-value = 0.99
## alternative hypothesis: two.sided
## 
## 
##   Outlier test:

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  sim
## outliers at both margin(s) = 1, observations = 47, p-value = 0.08972
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0005385317 0.1129377171
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                            0.0212766 
## 
## 
## ========== Diagnostics for NoWound model ==========

## 
##   Dispersion test:

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.97448, p-value = 0.934
## alternative hypothesis: two.sided
## 
## 
##   Outlier test:

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  sim
## outliers at both margin(s) = 0, observations = 47, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.00000000 0.07548573
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                                    0

cat("\n✓ All diagnostic plots saved to figures/diagnostics/\n")

## 
## ✓ All diagnostic plots saved to figures/diagnostics/

Predicted vs Observed

# Add predictions to data
percent_closure$pred_full <- predict(rate_full)

# Create scatterplot
p_fit <- ggplot(percent_closure, aes(x = pred_full, y = healing_rate)) +
  geom_point(color = "steelblue", size = 2, alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "darkred") +
  labs(
    x = "Predicted Healing Rate",
    y = "Observed Healing Rate",
    title = "Model Fit Check: rate_full"
  ) +
  theme_bw(base_size = 14)

print(p_fit)

ggsave(
  here("figures", "diagnostics", "fit_scatterplot.png"),
  p_fit,
  width = 10, height = 8, dpi = 300, bg = "white"
)

Residual Diagnostic Panel

library(patchwork)

# Residuals vs index
p_resid_plot <- ggplot(
  data.frame(resid = residuals(rate_full, type = "pearson")),
  aes(x = seq_along(resid), y = resid)
) +
  geom_point(alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  labs(title = "Pearson Residuals", x = "Observation", y = "Residual")

# Residuals vs fitted
p_fit_plot <- ggplot(
  data.frame(fitted = fitted(rate_full),
             resid = residuals(rate_full, type = "pearson")),
  aes(x = fitted, y = resid)
) +
  geom_point(alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  labs(title = "Residuals vs Fitted", x = "Fitted values", y = "Residuals")

# Combine panels
p_resid_panel <- (p_resid_plot | p_fit_plot) +
  plot_annotation(
    title = "Residual Diagnostics for Healing Rate Models",
    theme = theme(plot.title = element_text(size = 18, face = "bold"))
  )

print(p_resid_panel)

# Save panel
ggsave(
  here("figures","diagnostics","healing_rate_diagnostics_panel.png"),
  p_resid_panel,
  width = 14, height = 7, dpi = 300, bg = "white"
)

QQ-Plot for Normality

p_qq <- ggplot(
  data.frame(resid = residuals(rate_full, type = "pearson")),
  aes(sample = resid)
) +
  stat_qq(color = "steelblue", size = 2) +
  stat_qq_line(color = "darkred", linewidth = 1) +
  theme_bw(base_size = 14) +
  labs(
    title = "Normal Q-Q Plot of Residuals",
    x = "Theoretical Quantiles",
    y = "Sample Quantiles"
  )

print(p_qq)

ggsave(
  here("figures", "diagnostics", "qq_plot_residuals.png"),
  p_qq,
  width = 8, height = 8, dpi = 300, bg = "white"
)

Scale-Location Plot

p_scale_loc <- ggplot(
  data.frame(
    fitted = fitted(rate_full),
    sqrt_abs_resid = sqrt(abs(residuals(rate_full, type = "pearson")))
  ),
  aes(x = fitted, y = sqrt_abs_resid)
) +
  geom_point(alpha = 0.6, color = "steelblue", size = 2) +
  geom_smooth(se = TRUE, color = "darkred", method = "loess") +
  theme_bw(base_size = 14) +
  labs(
    title = "Scale-Location Plot",
    x = "Fitted Values",
    y = expression(sqrt("|Standardized Residuals|"))
  )

print(p_scale_loc)

ggsave(
  here("figures", "diagnostics", "scale_location_plot.png"),
  p_scale_loc,
  width = 10, height = 7, dpi = 300, bg = "white"
)

Diagnostic Summary Statistics

cat("\n=== Model Diagnostics Summary ===\n")

## 
## === Model Diagnostics Summary ===

cat("\nResidual standard deviation:", sigma(rate_full), "\n")

## 
## Residual standard deviation: 0.04185224

cat("Number of observations:", nobs(rate_full), "\n")

## Number of observations: 47

cat("\nResidual summary statistics:\n")

## 
## Residual summary statistics:

print(summary(residuals(rate_full, type = "pearson")))

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.114238 -0.022513  0.001759  0.000000  0.023531  0.089390

# Check for influential observations
cat("\n\n=== Checking for Influential Observations ===\n")

## 
## 
## === Checking for Influential Observations ===

influence_measures <- influence(rate_full, obs = TRUE)
cat("Cook's distance threshold (4/n):", 4/nobs(rate_full), "\n")

## Cook's distance threshold (4/n): 0.08510638

---

8. Session Information

sessionInfo()

## R version 4.5.1 (2025-06-13)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.0.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/Los_Angeles
## tzcode source: internal
## 
## attached base packages:
## [1] parallel  splines   stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] gamlss_5.5-0        nlme_3.1-168        gamlss.dist_6.1-1  
##  [4] gamlss.data_6.0-6   showtext_0.9-7      showtextdb_3.0     
##  [7] sysfonts_0.8.9      gt_1.0.0            patchwork_1.3.1    
## [10] viridis_0.6.5       viridisLite_0.4.2   ggpubr_0.6.1       
## [13] ggeffects_2.3.1     broom.mixed_0.2.9.6 parameters_0.28.0  
## [16] MuMIn_1.48.11       DHARMa_0.4.7        glmmTMB_1.1.12     
## [19] lme4_1.1-37         Matrix_1.7-3        janitor_2.2.1      
## [22] here_1.0.1          lubridate_1.9.4     forcats_1.0.0      
## [25] stringr_1.5.1       dplyr_1.1.4         purrr_1.1.0        
## [28] readr_2.1.5         tidyr_1.3.1         tibble_3.3.0       
## [31] ggplot2_3.5.2       tidyverse_2.0.0    
## 
## loaded via a namespace (and not attached):
##  [1] Rdpack_2.6.4        gridExtra_2.3       sandwich_3.1-1     
##  [4] rlang_1.1.6         magrittr_2.0.3      multcomp_1.4-28    
##  [7] snakecase_0.11.1    furrr_0.3.1         compiler_4.5.1     
## [10] mgcv_1.9-3          systemfonts_1.2.3   vctrs_0.6.5        
## [13] crayon_1.5.3        pkgconfig_2.0.3     fastmap_1.2.0      
## [16] backports_1.5.0     labeling_0.4.3      utf8_1.2.6         
## [19] rmarkdown_2.29      tzdb_0.5.0          nloptr_2.2.1       
## [22] ragg_1.4.0          bit_4.6.0           xfun_0.53          
## [25] cachem_1.1.0        jsonlite_2.0.0      broom_1.0.9        
## [28] gap.datasets_0.0.6  R6_2.6.1            bslib_0.9.0        
## [31] stringi_1.8.7       RColorBrewer_1.1-3  parallelly_1.45.1  
## [34] car_3.1-3           boot_1.3-31         jquerylib_0.1.4    
## [37] numDeriv_2016.8-1.1 estimability_1.5.1  Rcpp_1.1.0         
## [40] knitr_1.50          zoo_1.8-14          timechange_0.3.0   
## [43] tidyselect_1.2.1    abind_1.4-8         yaml_2.3.10        
## [46] TMB_1.9.17          codetools_0.2-20    curl_7.0.0         
## [49] listenv_0.9.1       lattice_0.22-7      withr_3.0.2        
## [52] bayestestR_0.16.1   coda_0.19-4.1       evaluate_1.0.4     
## [55] future_1.67.0       survival_3.8-3      xml2_1.4.0         
## [58] pillar_1.11.0       gap_1.6             carData_3.0-5      
## [61] stats4_4.5.1        reformulas_0.4.1    insight_1.4.0      
## [64] generics_0.1.4      vroom_1.6.5         rprojroot_2.1.0    
## [67] hms_1.1.3           scales_1.4.0        minqa_1.2.8        
## [70] globals_0.18.0      xtable_1.8-4        glue_1.8.0         
## [73] emmeans_1.11.2      tools_4.5.1         ggsignif_0.6.4     
## [76] fs_1.6.6            mvtnorm_1.3-3       grid_4.5.1         
## [79] rbibutils_2.3       datawizard_1.2.0    Formula_1.2-5      
## [82] cli_3.6.5           textshaping_1.0.1   gtable_0.3.6       
## [85] rstatix_0.7.2       sass_0.4.10         digest_0.6.37      
## [88] TH.data_1.1-3       farver_2.1.2        htmltools_0.5.8.1  
## [91] lifecycle_1.0.4     bit64_4.6.0-1       MASS_7.3-65

---

Analysis Complete ✓

All figures and tables saved to figures/ directory.